Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add corfidi MCS motion #3116

Merged
merged 12 commits into from
Nov 15, 2023
Merged

add corfidi MCS motion #3116

merged 12 commits into from
Nov 15, 2023

Conversation

wx4stg
Copy link
Contributor

@wx4stg wx4stg commented Jul 13, 2023

This adds support for Mesoscale Convective System storm motion as described in Corfidi 2003: https://journals.ametsoc.org/view/journals/wefo/18/6/1520-0434_2003_018_0997_cpampf_2_0_co_2.xml

Note that this implementation differs slightly from that of sharppy/SHARPpy's slightly, in that they use a non-pressure weighted wind for the cloud layer winds (this shouldn't have much of an effect on output in most cases) and they take the low-level jet as the mean of 0->1.5km AGL winds. My implementation looks for the maximum below 1.5 km AGL as I believe this is more in line with the original authors' intentions, "In the absence of a distinct low-level speed maximum in the vertical direction, the strongest wind in the lowest 5000 ft (1.5 km) is generally used, in accordance with Bonner (1968)". I'm extremely open to discussion on whether this is actually "better" or not though.

This is also my first contribution, so feel free to edit my documentation/code style/etc.

Thanks!

  • Tests added
  • Fully documented

@wx4stg wx4stg requested a review from a team as a code owner July 13, 2023 19:47
@wx4stg wx4stg requested review from dopplershift and removed request for a team July 13, 2023 19:47
@CLAassistant
Copy link

CLAassistant commented Jul 13, 2023

CLA assistant check
All committers have signed the CLA.

@wx4stg
Copy link
Contributor Author

wx4stg commented Jul 15, 2023

seems I didn't completely cover this in my unit test, stand by...

@wx4stg
Copy link
Contributor Author

wx4stg commented Jul 16, 2023

I'm on the fence as to whether I should be using pressure_weighted_mean or weighted_continuous average for the "cloud-layer winds". The original paper describes "a vector that represents cell advection by the mean cloud-layer wind (with “cloud layer” taken to be the 850–300-hPa layer)1" (Corfidi 2003) but then uses (wind_at_850+wind_at_700+wind_at_500+wind_at_300)/4 to estimate this.

Thinking more on this, this seems to weight the winds at higher pressures more, so maybe I should switch? It doesn't make a huge difference...

@dopplershift
Copy link
Member

@wx4stg Just wanted to let you know that I have been looking at this when I have the cycles and been trying to figure out how I want to proceed on the "cloud-layer winds" problem.

@wx4stg
Copy link
Contributor Author

wx4stg commented Aug 11, 2023

yeah--no problem, I meant to attend last week's dev call but forgot (oops), but I plan to start attending those once the semester gets going.

@dopplershift dopplershift added this to the September 2023 milestone Aug 30, 2023
@kgoebber
Copy link
Collaborator

Okay, so I've dug into this a bit and have attempted to reproduce the values produced in the paper and have not been able to do so (although very close in most circumstances), even with implementing their form of averaging for the cloud layer wind. My best guess is that the LLJ wind is the cause of the differences I'm seeing. I was not able to get the same value they propose for the 16 August 1997 case using either the 850-hPa LLJ method or the maximum wind in the lowest 1500 m.

My thoughts at this point are that the pressure weighted average is likely better. SharpPy averages the winds over the lower 1500 m for the LLJ and also implements non-pressure weighted averaging for both the LLJ and the cloud layer winds.

@mwilson14 Thoughts since you have thought about a lot of sounding indices in the past.

@kgoebber
Copy link
Collaborator

kgoebber commented Oct 6, 2023

Okay, here is the latest of what I have found when working to compare to the table provided in Confide (2003) with the current version of the implemented vector in the PR. In general, it doesn't match the values in the paper, but is close for some soundings, reasonable for many more, and closer to the observed values in other instances. I tried a couple of variations with only using the four levels and using a mean wind for the LLJ and all of those combinations made the errors greater. The only difference that was fairly small and zero in a number of cases was using the 850-hPa wind for the LLJ instead of the max wind in the lowest 1500 m. My best guess as to the differences is that unique LLJ vectors were used for different cases.

I think we should implement a method within the function to be able to specify a particular LLJ u and v-component to override the default behavior. I think it is an open conversation whether to use the maximum in the lowest 1500 m or just the 850-hPa wind for the LLJ. If we do the latter, then we could reduce the need for the height as an extra input variable - though this would not be fully in line with the what the paper describes as a best practice.

Note: I used Wyoming Sounding Archive and had to change the station name for

  • 87081912 OMA to OVN
  • 89070212 OKC to OUN
  • 89071712 OKC to OUN
  • 90052800 SIL to LCH
  • 91040912 LIT to LZK
Date Station Forecast wdir/spd Obs wdir/spd MetPy wdir/spd Diff w/ MetPy wdir Diff w/ MetPy spd FCST-Obs wdir FCST-Obs spd
83052000 VCT 235/60 260/55 229/58 -31 3 -25 5
83052000 DRT 245/45 265/40 248/57 -17 17 -20 5
83070112 PIA 270/43 290/45 242/34 -48 -11 -20 -2
83071912 BIS 290/57 300/50 264/62 -36 12 -10 7
83072000 GRB 300/54 310/50 296/59 -14 9 -10 4
83072200 WAL 315/53 315/45 308/72 -7 27 0 8
87061912 AMA 295/30 315/30 270/26 -45 -4 -20 0
87070412 UMN 280/30 290/30 282/30 -8 0 -10 0
87070512 UMN 270/32 270/35 306/26 36 -9 0 -3
87071112 STC 225/48 230/44 224/27 -6 -17 -5 4
87072012 STC 225/41 290/35 240/49 -50 14 -65 6
87081912 OVN 290/55 320/30 284/63 -36 33 -30 25
87082000 UMN 310/40 330/30 307/42 -23 12 -20 10
87111600 LCH 260/07 260/10 231/54 -29 44 0 -3
88040600 SLO 255/49 270/40 252/49 -18 9 -15 9
88051000 HTS 280/72 270/30 278/69 8 39 10 42
88071700 ACY 325/24 315/22 340/36 25 14 10 2
88081600 WAL 310/25 300/30 326/26 26 -4 10 -5
88081800 ACY 330/57 330/55 334/60 4 5 0 2
89022112 CHS 230/51 230/37 242/56 12 19 0 14
89050500 SEP 290/53 310/65 288/51 -22 -14 -20 -12
89050600 CHS 240/45 230/45 258/50 28 5 10 0
89061700 WAL 215/46 220/45 225/45 5 -0 -5 1
89070212 OUN 340/50 340/50 337/54 -3 4 0 0
89071712 OUN 300/32 315/35 300/26 -15 -9 -15 -3
89080600 HTS 270/30 290/25 299/18 9 -7 -20 5
89080700 AMA 330/25 350/25 326/29 -24 4 -20 0
89112100 IAD 285/60 285/55 289/86 4 31 0 5
90021012 AHN 240/57 270/50 226/82 -44 32 -30 7
90052800 LCH 250/30 270/30 277/67 7 37 -20 0
90062000 TOP 290/40 270/45 281/45 11 0 20 -5
90082612 GRB 260/27 300/30 241/31 -59 1 -40 -3
91040912 LZK 240/47 240/40 249/42 9 2 0 7
91050500 GGG 240/38 250/30 234/41 -16 11 -10 8
91050600 AYS 270/35 270/40 272/33 2 -7 0 -5
91050700 ACY 240/31 270/40 230/59 -40 19 -30 -9
91052900 FNT 290/37 300/35 296/39 -4 4 -10 2
91070800 FNT 260/50 270/45 258/49 -12 4 -10 5
92070300 FNT 300/48 290/40 301/43 11 3 10 8
92070300 UMN 250/47 270/45 241/43 -29 -2 -20 2
93060412 PAH 270/50 280/52 256/41 -24 -11 -10 -2
93060500 HAT 285/67 300/60 274/67 -26 7 -15 7
93080100 PAH 330/38 320/40 322/25 2 -15 10 -2
93080100 TOP 315/35 320/35 305/44 -15 9 -5 0
94041800 PIA 300/55 280/55 312/69 32 14 20 0
96050512 LZK 270/40 270/70 252/34 -18 -36 0 -30
96052112 CHH 270/50 270/60 263/48 -7 -12 0 -10
97081612 DTX 265/40 275/45 255/51 -20 6 -10 -5

@mwilson14
Copy link
Contributor

This implementation makes sense to me! I agree that the pressure-weighted wind is likely better and that having an option in the function to specify a LLJ vector is a good idea. Also, sorry for missing that first comment I was tagged in, that was the weekend I moved to Boulder so things were a bit chaotic.

@kgoebber
Copy link
Collaborator

@wx4stg, okay so I think we have a final recommendation on implementation for this PR!

  1. Add in an input option for a u and v LLJ. This will allow a more nimble method for end users who might know or want to specify a certain LLJ instead of using the default option that is implemented in the function.
  2. Instead of using the maximum in the lowest 1500 m, to use the maximum wind from the surface to 850 hPa. This would allow us to reduce the number of input variables for users and not change the default value coming out of the function too much.
  3. Put a second example in the function that demonstrates finding the maximum LLJ in the lowest 1500 m and putting that LLJ value in the optional variables.

Let @dcamron or I know if you would like any help or want to punt any of these updates to us if you don't have the cycles at the midpoint of the current semester.

@wx4stg
Copy link
Contributor Author

wx4stg commented Oct 19, 2023

@kgoebber thanks for the review! I have a few questions..

Add in an input option for a u and v LLJ. This will allow a more nimble method for end users who might know or want to specify a certain LLJ instead of using the default option that is implemented in the function.

Ok, just to clarify, I'll add keyword args, something like llj_u and llj_v, which are None by default. If the function is called without specifying, then proceed to find LLJ as maximum below 850 hPa?

Instead of using the maximum in the lowest 1500 m, to use the maximum wind from the surface to 850 hPa. This would allow us to reduce the number of input variables for users and not change the default value coming out of the function too much.

I agree with the idea of using pressure coordinate height for this, but this introduces a bit of an edge case: what if the sounding is from something like Denver where the surface pressure occasionally dips below 850 hPa? Raise ValueError? High altitude MCSes are... certainly not very typical, but just curious about how to handle this.

Let dcamron or I know if you would like any help or want to punt any of these updates to us if you don't have the cycles at the midpoint of the current semester.

Hilariously, these updates could not have come at a better time, we just finished midterms so I have a little breathing room this week :)

@kgoebber
Copy link
Collaborator

@kgoebber thanks for the review! I have a few questions..

Add in an input option for a u and v LLJ. This will allow a more nimble method for end users who might know or want to specify a certain LLJ instead of using the default option that is implemented in the function.

Ok, just to clarify, I'll add keyword args, something like llj_u and llj_v, which are None by default. If the function is called without specifying, then proceed to find LLJ as maximum below 850 hPa?

Yep! So in the end it should be something like...

def corfidi_storm_motion(pressure, u, v, *, u_llj, v_llj):

@dopplershift, @dcamron thoughts on naming of LLJ u and v components?

Instead of using the maximum in the lowest 1500 m, to use the maximum wind from the surface to 850 hPa. This would allow us to reduce the number of input variables for users and not change the default value coming out of the function too much.

I agree with the idea of using pressure coordinate height for this, but this introduces a bit of an edge case: what if the sounding is from something like Denver where the surface pressure occasionally dips below 850 hPa? Raise ValueError? High altitude MCSes are... certainly not very typical, but just curious about how to handle this.

I think raise a value error would be good and prompt to specify an LLJ u and v component. @dopplershift @dcamron thoughts?

Let dcamron or I know if you would like any help or want to punt any of these updates to us if you don't have the cycles at the midpoint of the current semester.

Hilariously, these updates could not have come at a better time, we just finished midterms so I have a little breathing room this week :)

Excellent! Hope midterms went well!

@dopplershift
Copy link
Member

def corfidi_storm_motion(pressure, u, v, *, u_llj, v_llj):

@dopplershift, @dcamron thoughts on naming of LLJ u and v components?

While we usually try to avoid terse/abbreviated names so that the function is more self-documenting, u_low_level_jet seems a bit verbose and ugly for something we are requiring as a keyword argument, so I think I'm ok with eg. u_llj--we just of course have to spell out "low-level jet" when documenting what these are.

Instead of using the maximum in the lowest 1500 m, to use the maximum wind from the surface to 850 hPa. This would allow us to reduce the number of input variables for users and not change the default value coming out of the function too much.

I think raise a value error would be good and prompt to specify an LLJ u and v component. @dopplershift @dcamron thoughts?

Agreed. Something like ValueError("Unable to automatically determine low-level jet winds, please pass them in manually with u_llj and v_llj.)

@wx4stg
Copy link
Contributor Author

wx4stg commented Oct 30, 2023

@kgoebber I think this covers the requested changes and provides tests, if you still have the code from your earlier comparison, you can try it again.

strangely I can't get the bots to run..

@kgoebber
Copy link
Collaborator

There is a conflict with the main branch. There are a couple of ways to resolve the conflict. In general it is best of you update your main MetPy branch

git checkout main

git pull upstream main

git checkout <branch>

git rebase -i main

The last command is the rebase in an interactive mode. It will notify you if there are conflicts that are not resolved automatically that might need your intervention. Once you have rebased, then you can push the branch updates as per usual.

@wx4stg
Copy link
Contributor Author

wx4stg commented Oct 31, 2023

I see. I wasn't sure if y'all would want to handle that or not. thanks.

@dopplershift dopplershift merged commit f15ec65 into Unidata:main Nov 15, 2023
34 checks passed
@dopplershift dopplershift added Type: Feature New functionality Area: Calc Pertains to calculations labels Nov 15, 2023
@dopplershift
Copy link
Member

Thanks for working with us on getting this in @wx4stg !

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Area: Calc Pertains to calculations Type: Feature New functionality
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants